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Many difficult computational problems involve the simultaneous satisfaction of multiple con- 
straints which are individually easy to satisfy. Such problems occur in diffractive imaging, protein 
folding, constrained optimization (e.g., spin glasses), and satisfiability testing. We present a simple 
geometric framework to express and solve such problems and apply it to two benchmarks. In the first 
application (3SAT, a boolean satisfaction problem), the resulting method exhibits similar perfor- 
mance scaling as a leading context-specific algorithm (walksat). In the second application (sphere 
packing), the method allowed us to find improved solutions to some old and well-studied optimiza- 
tion problems. Based upon its simplicity and observed efficiency, we argue that this framework 
provides a competitive alternative to stochastic methods such as simulated annealing. 



Difficult problems can often be broken down into a col- 
lection of smaller, more tractable, subproblems. This is 
the basis of the divide and conquer approach, which ap- 
plies when the initial problem and the subproblems have 
a similar structure, and the global solution can be re- 
trieved from the solutions to the subproblems. Divide 
and conquer as a rule leads to very efficient algorithms. 
However, many difficult problems do not fit such an effi- 
cient framework. 

For example, consider the problem of determining the 
three-dimensional structure of a complex molecule given 
clues about the distances between particular pairs of 
atoms (from knowledge of chemical bonds, NMR mea- 
surement, etc.). As subproblems we might consider the 
substructures formed by small groups of atoms, since 
finding substructures satisfying local constraints is usu- 
ally not challenging. However, the location and orienta- 
tion in space of the substructures depends intricately and 
sensitively on their collective arrangement. Because the 
division into subproblems in this case does not lead to 
a practical algorithm, molecular geometry problems are 
usually transformed into optimization problems through 
the definition of a global cost function, and are then 
solved through stochastic optimization methods such as 
simulated annealing. 

In this Letter we introduce a general method for solv- 
ing constraint problems that takes advantage of the di- 
vision into subproblems. In broad terms the method dif- 
fers from stochastic searches in that the configurations 
explored are generated iteratively and deterministically. 
Each iterative step is defined by two fundamental opera- 
tions. In the first operation, the problem is divided into 
its constituent constraints, which are then solved inde- 
pendently, ignoring possible conflicts between different 
constraints. In the second operation, conflicts between 
constraints are resolved regardless of the satisfaction of 
the constraints. By a judicious application of these two 
operations, we obtain a search strategy which, at each 
step, solves all the subproblems separately and at the 
same time seeks to resolve conflicts between their solu- 
tions. We call this method divide and concur (D&C). 



The D&C approach can be applied to a wide range of 
problems, both discrete and continuous. We first show 
how D&C is applied to the boolean satisfiability problem 
(SAT), a standard benchmark in computer science. In 
this problem, the D&C approach exhibits similar scaling 
behavior to walksat, a leading SAT solver which out- 
performs general-purpose algorithms such as simulated 
annealing !1 . As a second example we study continu- 
ous sphere packing problems, which are formally simi- 
lar to the molecular geometry example mentioned above. 
The D&C approach matched or improved upon the best 
known packings in some well-studied, two-dimensional 
problems. In 10 dimensions it also discovered an inter- 
esting new sphere arrangement related to quasicrystals. 
The D&C approach therefore combines the advantages of 
general purpose algorithms (versatility, simplicity) with 
the performance of special purpose algorithms (such as 
walksat) . 

In D&C the individual constraints are first expressed 
as subsets of a Euclidean space K, thereby transforming 
the constraint satisfaction problem into the geometrical 
problem of finding a point in the intersection of multiple 
sets. The Euclidean space provides the setting to define 
distance-minimizing projections to each of the N con- 
straint sets. The projection operators {Pj}i = i...jv wm be 
the building blocks of the algorithm. Starting from an 
initial guess, one uses the projections to probe the con- 
straint sets and update the guess. This idea has been 
studied extensively in the context of convex constraint 
sets 0. 

Given N primary constraints expressed as subsets of 
K, we first define the product space K N , consisting of 
N copies (or replicas) of K |15j . We then define, in 
the product space, the 'divide' constraint D (enforcing 
one primary constraint on each replica) and the 'concur' 
constraint C (enforcing replica concurrence) [H [3 . The 
associated projections, acting ony = ® x^ 2 ) ® ■ ■ ■ <£> 
x^, are 

P D (y) = P x (x«) ® P 2 (x( 2 >) ® •• • ® Pjv(xW), (1) 
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which acts separately on each of the replicas, and 

Pc{y) = x®x®---®x, (2) 

which replaces the value of each replica by the average 
value x of all the replicas. In defining the concurrence 
projection, different weights Aj may be assigned to differ- 
ent constraints, i.e., x = Si(^ x ^)/ Si ^« [IS]- Chang- 
ing the weights is equivalent to changing the metric of 
the product space; this possibility will prove beneficial 
even in problems where all the constraints are formally 
equivalent. Note that both projections, Pc and Pd, act 
in a highly parallel sense: either by treating indepen- 
dently each replica (Pd), or by treating independently 
each variable across the replicas (Pc)- 

Through the product space construction the original 
constraint problem has been expressed as the problem of 
finding a point in the intersection of two sets, both of 
which have easily implemented projection operators. To 
proceed, we need a search strategy that can use a pair of 
projection operators (P a and Pb) to seek the intersection 
of two sets. The simplest approach is the alternating pro- 
jections scheme, where y„+i = P<j(P&(y n )) [HE]- Despite 
its success with convex constraints and some nonconvex 
problems, the alternating projections scheme is prone to 
getting stuck at fixed points which do not correspond to 
solutions. 

The difference map (DM) is an improvement upon al- 
ternating projections which emerged in response to the 
nonconvex constraints arising in diffractive imaging (the 
phase problem) [51 [BJ. It is defined by a slightly more 
elaborate set of rules, namely 

y„+i =y„ + (3 (P a o f b (y n ) -P b o f a (y n )) 
fi(y n ) =(i + it)Pi(yn) - liy-n i = a, b, 

with 7 a = —1//3 and jb = 1/P- The parameter f3 can 
have either sign and is chosen to improve performance. 
If the iteration reaches a fixed point y*, an intersection 
of the constraint sets has been found. The solution y so i is 
obtained from the fixed point using y so i — P a ° fb(y*) — 
Pb ° fa(y*)- Note that the fixed point itself is not neces- 
sarily a solution; in fact, there typically is a continuum 
of fixed points associated with every solution. 

The DM was recently used to solve a variety of difficult 
computational problems, including protein folding [7, 8J, 
boolean satisfiability, Diophantinc equations, graph col- 
oring, and spin glasses [§]• Most of these applications 
relied on a decoupling of constraints through the use of a 
dual set of variables, as in linear programming. The di- 
vide and concur approach we introduce here, defined by 
the use of the difference map with P a — Pc and Pb = Pd , 
is significantly more versatile and systematic. 

The boolean satisfiability problem 3SAT is one of the 
most extensively studied problems in constraint satisfac- 
tion. The challenge is to find an assignement for N v 



boolean variables that satisfies a list of N c boolean con- 
straints, or clauses. Each clause is an OR statement in- 
volving three literals, i\ V £2 V £3, where each literal li 
represents either one of the N v boolean variables or its 
negation. 

A D&C formulation of 3SAT is obtained by associating 
a real-valued search variable to each 3SAT literal, where 
the values {1,-1} are taken to mean {True, False}. The 
constraint D requires that each clause is satisfied; that is, 
each literal must have value ±1, with at least one literal 
per clause having value 1. In other words, to each clause 
corresponds a variable triplet, which is projected by Pd 
to the nearest of the seven satisfying assignments for this 
clause. Geometrically, these correspond to seven vertices 
of a cube. In this application Pc ensures that all literals 
associated with the same boolean variable concur (with 
due regard to negations). Since each constraint (clause) 
involves only three variables, the reduced search space[15] 
has dimension 3iV c . For simplicity, we give equal weight 
to each constraint (A, = 1). 

We compared the performance of the D&C algorithm 
with WALKSAT [T] on a collection of 3SAT problem in- 
stances ranging from N v = 50 to N v = 25600, with fixed 
ratio a = N c /N v = 4.2, a value for which randomly gen- 
erated instances are expected to be difficult [10] . Random 
instances were generated using the program makewff 
(distributed with walksat), and instances that were not 
solved by either walksat or the D&C algorithm were 
discarded. Each algorithm was applied 10 times to each 
instance, starting from different random initial condi- 
tions. The median number of variable updates required 
to find the solution is plotted in Figure [T] The number 
of variable updates in walksat equals the total number 
of flips of the boolean variables. In the D&C algorithm 
it is the total number of nonzero updates of any of the 
real- valued search variables (literals). 

Fig. [l] shows that Walksat (with the 'noise' parame- 
ter fixed at the value p = 0.57) and the D&C algorithm 
(with (3 — 0.9) have similar performance behavior. Not 
only do they both find the same problems easy and the 
same problems hard (which is not unexpected), but the 
scaling of the number of variable updates needed to reach 
the solution, as a function of problem size, is also similar. 
Such a similarity is surprising, considering the difference 
in search strategies. 

walksat uses pseudorandom processes (or 'noise') to 
update the variables asynchronously. In D&C, on the 
other hand, the update rule is completely deterministic 
and is applied synchronously to many variables. Fig. 
1 also shows that choosing suboptimal parameters for 
either algorithm results in rapid performance degrada- 
tion for large problem sizes. Even though the scaling of 
the variable updates are similar for walksat and D&C, 
our implementation of D&C required significantly more 
CPU time (between 4 and 200 times, depending on the 
instance) than walksat. Work on an optimized im- 
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FIG. 1: Median number of variable updates needed to find a 
solution for walksat (WS) and divide and concur (D&C) on 
the same set of random 3SAT instances with a — 4.2. Each 
median was calculated by solving the same instance 10 times 
starting from different random initial guesses, for parameter 
values P = 0.9 (D&C) and p = 0.57 (walksat). Variations 
resulting from changing (3 and p are indicated by the shaded 
areas; both methods exhibit parameter sensitivity for prob- 
lems with more than 10 4 variables. A point at the top edge 
indicates that the median exceeded the cutoff on the number 
of updates, 3 x 10 10 . 

plementation of the D&C algorithm is in progress and 
should allow an easier exploration of the behavior of the 
method for larger problem sizes. 

Another constraint problem which has been exten- 
sively studied is the packing of n spheres in a finite D- 
dimensional volume (see, e.g., Refs. jTTJ [TSJ [T3] and 
references therein). The constraint formulation of this 
problem is more directly geometrical than boolean sat- 
isfiability. Since each sphere must avoid n — 1 other 
spheres and lie within a certain volume, there are al- 
together n constraints per sphere. The reduced search 
space [15] requires one _D-dimensional variable replica for 
every sphere participating in a constraint, for a net search 
space dimensionality of Dn 2 . 

Within the framework of D&C there is a formal simi- 
larity in the constraint structure of packing spheres and 
3SAT. Just as every boolean variable is constrained by 
each of the clauses where it occurs, every sphere in a 
packing has a volume exclusion relationship with each of 
the other spheres in the packing: ||x a — x/,| > m a b- This 
similarity and the success of D&C with 3SAT is strong 
motivation to apply D&C to the sphere packing problem. 

Near the solution of any n-sphere packing problem, 
the number of relevant exclusion constraints (contact- 
ing pairs) grows only as n (for fixed D) while the to- 
tal number of constraints is 0(n 2 ). In the D&C ap- 
proach it is possible to increase the weight of these rel- 
evant pairs by dynamically adjusting the corresponding 
metric weight A a ;,. This results in considerable perfor- 
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FIG. 2: An example of an improved packing for 169 disks in 
a square found by the D&C algorithm. The figure on the left 
shows the previously best known packing [T3], with density 
0.8393. The density of the improved packing shown on the 
right is 0.8399. Contacts are shown with dotted lines; colors 
indicate the number of contacts. 

mance improvement. At the end of each DM step we 
used Xab — > cA afc + (1 - a)exp(-adab), where d ab is 
the current distance between the pair[T7]. We used the 
value a = 0.99 to ensure that the metric update is quasi- 
adiabatic (i.e., slow on the time scale of variable updates), 
and a ~ 30. 

We first consider packings of n equal disks of diameter 
m in a unit square, and take as a starting point the best 
known packing diameters m* from Ref. |13j . This prob- 
lem is quite challenging, due to the coexistence of many 
different arrangements with similar density. We tested 
the D&C algorithm for each value of n in the range 2- 
200. For each n, we generated up to 400 random initial 
guesses. For each initial guess, a small value of the di- 
ameter m was chosen, and a packing was sought. When 
a solution was found, m was increased, and the process 
was repeated until the algorithm failed to find a pack- 
ing, or until the best known packing diameter m* was 
reached. In the latter case the target was increased be- 
yond m* with the hope of finding a denser packing. No 
information about the known packings was used, apart 
from their densities. 

For 143 of the 197 values of n a packing with diameter 
close to the optimal packing (m > m* — 10 -9 ) was found. 
More surprisingly, improved packings were found in 38 
cases. The smallest n for which an improved packing was 
found is 91. The largest improvement was for n = 182, 
for which a packing was found with m = m* + 4.6 x 10~ 5 . 
For 28 values of n a packing was found with m > m* + 
1 x 10~ 6 . An example of such an improved packing is 
shown in Figure [2] 

When packing many disks the optimization challenge is 
easy to identify as a contest between close-packing in the 
bulk and an efficient match to the boundary. In higher 
dimensions the structure of the solution is not so easily 
characterized, and we can look to the D&C method as 
an unbiased tool for exploration. A classic problem in 
geometry is to determine kissing numbers td- the max- 
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imum number of unit spheres that can be packed in D- 
dimensions, so that each contacts a given unit sphere. 
Early investigations of this problem were stimulated by 
a debate between Newton and Gregory, who disputed the 
value of T3. The only known kissing numbers are ti = 2, 
t-2 = 6, r 3 = 12, r 4 = 24, r 8 = 240, and r 24 = 196 5 60. 
In dimension 1-8, and also 16-24, the best known lower 
bounds on td are given by the number of minimal vectors 
in the unique laminated lattice of the same dimension[TT] . 
For dimension 9-15 the best bounds are obtained from 
constructions based on error-correcting codes [IT . Dis- 
coveries of novel packings in higher dimensions has for 
the most part been achieved through mathematical in- 
spiration. Unbiased searches, defined only by the basic 
constraints, have to our knowledge not been attempted 
beyond dimension 5[14 . This raises the possibility that 
interesting packings in high dimensions may have escaped 
detection only for lack of imagination. 

With minimal adjustment to the above procedure for 
finding disk packings, we were able to find kissing ar- 
rangements as good as the best known in dimension 2- 
4, 6, and 8. After introducing just the assumption of 
inversion symmetry, optimal packings were obtained in 
all dimensions up to 8. Our searches in higher dimen- 
sions have so far revealed an interesting new packing in 
dimension 10. It is easy to understand why this pack- 
ing was missed. Constructions based on integral lattices 
and error-correcting codes all have the property that the 
cosine of the angle subtended by any two spheres is ratio- 
nal. The packing of 378 spheres discovered by the D&C 
algorithm has all cosines in a set that includes irrational 
numbers: {±l,±l/2, (±3 ± %/3)/12,0}. An analysis of 
the coordinates obtained by the algorithm has revealed 
that these 378 sphere positions are expressible as unique 
integer multiples of a basis of 12 vectors. The construc- 
tion has a strong relationship to quasicrystals, where the 
excess dimension of the basis accounts for irrational rela- 
tionships in the geometry. The algorithm, of course, had 
no knowledge of quasicrystal geometry. 

This 'irrational' structure emerged as soon as the num- 
ber of spheres was increased above 372, the largest known 
kissing number for 10-dimensional lattices jTT]. The 
same irrational arrangement was also found for up to 
384 spheres; the 6 additional spheres were accomodated 
in holes of the structure (and have continuously variable 
cosines). Finally, the algorithm has so far been unsuc- 
cessful in discovering the best known kissing arrangement 
in 10 dimensions, with kissing number 500. 

The divide and concur approach provides a natural 
framework in which to address various hard computa- 
tional problems. In two benchmark applications, 3SAT 
and sphere packing, the D&C approach compares with, 
and in some cases improves upon, state-of-the-art spe- 
cialized methods. The uniform mechanism provided by 
the difference map, for finding solutions given a pair of 
constraint projections, makes the D&C approach almost 



as easy to implement as general-purpose sampling algo- 
rithms such as simulated annealing. Most of the problem- 
specific development needed, in this framework, is the 
definition of the appropriate projection operators. We 
believe the latter are able to exploit important elements 
of the problem structure not accessed by stochastic sam- 
pling, and that this accounts for the superior performance 
of D&C. 
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